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PREFACE 
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responsible  for  publishing  revisions  or  for  identifying  the  obsoles- 
cence of  its  technical  publications. 

PURPOSE:  To  contribute  technical  information  to  the  field  of  geodesy 

by  describing  a method  for  estimating  gravity  anomalies  from  satellite 
derived  geoid  height  data. 
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DISTRIBUTION:  Qualified  requesters  may  obtain  copies  from  the  Defense 
Documentation  Center,  Cameron  Station,  Virginia  22314. 
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ABSTRACT 

A gridded  data  base  of  geoid  heights  can  be  developed  from 
GEOS-III  satellite  altimetry  data.  The  derivation  of  two  methods  is 
discussed.  A procedure  for  estimating  gravity  anomalies  from  the 
resulting  geoid  heights  is  described  with  simulations  shown  in 
tabular  and  graphic  form.  The  report  also  discusses  the  computation 
of  the  calibration  geoid  for  the  GEOS-III  test  region,  the  data  base 
from  altimetry  data  and  the  simulation  of  test  data  for  the  estima- 
tion program.  The  method  of  localization  for  estimation  and  the 
methods  of  testing  the  computer  program  are  given  in  some  detail. 
Test  results  which  were  computed  from  mld-1974  to  early  1975  are 
Included  along  with  supporting  tables  and  figures. 
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I.  INTRODUCTION 

u Prior  to  the  launching  of  the  first  artificial  satellite  in 

1957,  the  best  approximation  for  the  gravitational  potential  field  of 
the  earth  was  a spherical  harmonic  series  to  degree  and  order  four 
obtained  from  studies  of  the  moon's  orbit  and  gravity  observations  at 
a few  scattered  land  stations.  Extensive  satellite  tracking  data  has 
made  it  possible  to  analyze  orbital  perturbations  for  a more  detailed 
■ representation  of  the  potential  field.  In  recent  years  technical 

developments  in  improved  satellite  tracking  and  advances  in  missile 
weapon  systems  have  increased  the  need  for  a more  accurate  represen- 
tation of  the  potential  field  to  higher  harmonic  degree  and  order. 

The  incorporation  of  orbital  data  with  an  increasing  number  of  gravity 
observations  on  land  has  allowed  approximation  of  the  potential  field 
to  degree  and  order  twenty  with  several  zonal  and  tesseral  coeffi- 
cients beyond  that  level.  The  models  developed  by  several  investiga- 
tors during  the  last  few  years  are  fairly  consistent  to  degree  and 
order  twelve.  Usually,  the  lower  degree  harmonics  are  based  on  orbital 
data  and  the  higher  degree  harmonics  are  developed  from  surface 
observations.  Development  of  harmonic  series  to  higher  degree  and 
order  requires  more  detailed  knowledge  of  the  gravity  field  over  the 
entire  surface  of  the  earth.  Gravity  data,  particularly  ocean  gravity 
data,  is  generally  not  adequate  outside  North  America. 

The  National  Aeronautics  and  Space  Administration  (NASA) 
included  a radar  altimeter  in  the  Geodynamics  Experimental  Ocean 
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Satellite  (GEOS- III),  launched  9 April  1975.  The  purpose  of  the  altim- 
eter is  to  measure  the  ocean  surface  topography  in  order  to  obtain  a 
detailed  representation  of  the  geoid  for  the  ocean  regions.  Geoid 
heights,  deflections  of  the  vertical,  gravity  anomalies  and  the  gravi- 
tational potential  are  related  in  such  a way  that  observations  which 
are  adequate  for  computing  mean  geoid  heights  are  also  sufficient 
for  computing  mean  values  of  the  other  geodetic  parameters.  Principal 
investigators  were  selected  by  NASA  to  analyze  and  evaluate  the  data 
collected  from  the  various  experiments  on  board  the  satellite. 

As  a user  of  gravity  information,  DoD  has  developed  a 
separate  exploitation  plan  [1]  to  utilize  GEOS- 1 1 1 altimetry  data  for 
geoid,  deflection  of  the  vertical  and  mean  gravity  anomaly  information. 
As  a part  of  this  plan,  the  Naval  Research  Laboratory  (NRL)  corrects 
the  altimetry  data  for  tides,  winds  and  ocean  currents.  The  Naval 
Surface  Weapons  Center,  Dahlgren  Laboratory  (NSWC  DL)  computes  the 
corrected  heights  of  the  satellite  above  the  geoid,  along  track 
geoid  heights,  geocentric  radius  of  the  satellite  position  and  the 
ellipsoidal  radius  of  the  suborbital  point  (Figure  1). 

This  data  is  then  provided  to  the  Defense  Mapping  Agency 
Aerospace  Center  (DMAAC)  for  development  of  a gridded  data  base  of 
estimated  1°  x 1°  and  5°  x 5°  mean  geoid  heights  and  gravity  anomalies. 
This  data,  when  combined  with  mean  anomalies  from  surface  gravity  data 
and  satellite  orbital  analysis,  will  be  used  to  derive  an  improved 
earth  gravitational  model. 
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The  purpose  of  this  report  is  to  show  that  estimation  of 
gravity  anomalies  from  gridded  geoid  height  data  can  be  effectively 


. • • 

accomplished  with  the  methods  described.  For  specified  accuracy 
levels  of  the  geoid  height  data,  an  indication  of  the  accuracy 

attainable  for  estimating  mean  gravity  anomalies  over  1°  x 1°  • • 

squares  with  an  unweighted  least  squares  estimator  is  presented. 
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II.  GRAVIMETRIC  GEOID  HEIGHTS 
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, To  provide  a standard  for  comparison  of  satellite  derived 

geoid  heights,  DMAAC  computed  gravimetric  geoid  heights  over  a region 
of  the  North  Atlantic  Ocean  designated  for  GEOS- III  calibration.  The 
calculation  required  a worldwide  integration  of  5‘  x 5',  15'  x 15', 
lu  x 1°  and  5°  x 5°  elements,  (Figure  2). 

The  mean  free  air  gravity  anomalies  and  their  accuracies 
for  the  test  region  are  stored  on  magnetic  tapes  in  the  DoD  Gravity 
Library.  These  tapes  are  continually  updated.  Surface  point  values 
are  not  available  for  all  the  required  5'  x 5'  elements.  Therefore, 
interim  values  have  been  obtained  for  the  remaining  elements  by  least 
squares  prediction.  The  interim  values  constitute  about  50  percent 
of  all  the  required  mean  5'  x 5'  gravity  anomalies  in  the  calibration 
[ area.  The  15'  x 15'  and  1°  x 1°  values  in  the  test  region  are  means 

of  5'  x 5'  elements.  The  values  for  the  worldwide  array  of  5°  x 5° 
elements  are  averages  of  1°  x 1°  elements  obtained  from  observed 
data,  geophysical  correlation,  extrapolation,  etc. 

Values  for  all  5°  x 5°  elements  are  referenced  to  the  World 
Geodetic  System  1972  (WGS  72)  [2].  The  anomaly  data  referenced  to 
the  WGS  72  gravity  formula  was  used  to  compute  the  geoid  heights  at 
the  centers  of  the  1°  x 1°  surface  elements  in  the  calibration  region. 
A contour  plot  of  these  geoid  heights  is  shown  in  Figure  3.  Computa- 
tions are  being  performed  to  expand  the  calibration  region  to  include 
all  of  the  North  Atlantic  to  65°N,  the  limit  of  GEOS-III  coverage. 
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Fioura  3.  Gravimetric  Geoid  Heights  for  GEOS-III 
Calibration  Test  Renion  (meters). 
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III.  THE  6RIDDED  DATA  BASE 

The  gridded  data  base  consists  of  geold  heights  and  corre- 
sponding error  variances  estimated  at  the  centers  of  1°  x 1°  elements 
from  GEOS- III  data  (Figure  4).  The  geoid  height  observations  are 
regarded  as  the  sum  of  a stochastic  signal  (true  geoid  height)  and  a 
stochastic  noise  (error).  This  section  discusses  two  methods  for 
obtaining  the  gridded  geoid  height  data  base:  point  prediction  and 

profile  prediction. 

A.  Point  Prediction  Method 

The  centered  1°  x 1°  values  in  the  gridded  data  base 
are  estimated  by  Wiener  filtering  and  prediction  which  is  a discrete 
least  squares  estimator  [3].  It  is  assumed  that  the  estimator  Np 
at  a point  P is  a linear  function  of  a set  of  observations  N..  with 
prediction  coefficients  ap.  such  that 

n 

Np  = <*pi  N1  + ap2  N2  + ...  aRn  Nn  = ^ api  ^ . (1) 

The  prediction  coefficients  api-  depend  on  the  geoid 
height  autocovariance  function  and  the  positions  of  the  n observation 
points  relative  to  the  calculation  point  P.  The  best  estimates  of 
the  coefficients  (from  observed  data)  are  obtained  by  minimizing  the 
mean  square  error  of  prediction.  If  Np  is  the  true  value  of  N at 
point  P and  ep  is  the  error  of  prediction,  then 

£p  = Np  - Np  = Np  - S api  Nr  (2) 


f 2P  = ^NP  " i^1  “pi  Ni^Np  " “Pi  V’ 


(3) 


and 


n n 


-2P  = N2p  - 2 ^ “Pi  Np  Ni  + “Pi  “Pk  NiNk‘ 


(4) 


The  geoid  height  covariances  are  estimated  by  averaging 
over  the  region  of  interest, 

M{N.N. } = C..  = the  covariance  between  the  observed  geoid  heights 


i k ik 


at  points  i and  k. 


M{NnN. } = Cn.  = the  covariance  between  the  geoid  heights  at  the 


Pi  Pi 


computation  point  P and  the  observed  point  i. 


M(N2p}  = CQ  = the  variance  of  the  geoid  height. 

The  mean  square  error  of  prediction  of  the  geoid  height  at  point  P is 

n n n 


M{e2p}  = m2p  = CQ  - 2 ^ ap.  Cp.  + ^ ^ aPi  apk  Cik* 


(5) 


The  mean  square  prediction  error  is  minimized  by  setting  the  deriva- 
tive with  respect  to  api-  to  zero  and  solving  for  the  coefficients. 
Thus, 


8nrD  n 

2 cpi  + 2 pk  c<k*° 


(1=1 . 2 ...  n),  (6) 


which  reduces  to  a system  of  n linear  equations 

n 

1 

k= 


^ Cik  aPk  = CPi * 


(7) 


with  apk  as  unknown  coefficients.  The  solution  for  the  coefficients 


is 


“Pk  = Cik_1  CPi ' 


(8) 
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Then  the  predicted  geoid  height  becomes 


• • 


n n 


NP  = ^ aPk  Nk  = 1J]  Cik_1  CPi  V 


which  in  matrix  form  is 


C11  C12  Cln 
C21  C22  •••  C2n 


-1  r n. 


• • 


Np  = (Cpl>  Cp2  ...  CRn) 


, 00) 


c , c „ c i 

nl  n2  nn 


where  N2>  ...,  Nn  are  the  observed  geoid  heights. 

The  error  variance  of  the  least  squares  predictor  is 
calculated  by  inserting  the  coefficients  (8)  into  the  error  equation 


(5)  which  becomes 


n n 


2P  = °2P  C0  "2  Cik_1  CPi 


+ i k j ? Cik'1  Cpi  Cpj  Ch‘  (11) 

with  the  subscripts  j and  I representing  the  same  set  of  observed 
values  as  i and  k. 


Since 


Vr  - 1 c =/l  i^j~k» 
\ Lji  hi  ' *0  if  j7k. 


• • 


the  error  equation  becomes 


m2p  = C0  ’ iJ1  Cik_1  CPi  Sk- 


in matrix  form  this  is 


C11  C12  • 
C21  C22  ' 


• Cln|  _1  |CPV 


m2p  * c0  " (CP1*  CP2  V 


[Si  Cn2  CnnJ  [_CpJ 

The  matrix  [C.. of  the  autocovariances  of  geoid  heights  at  the 
observed  points  is  used  to  predict  the  geoid  height  and  its  error 
variance.  The  estimators  (10)  and  (13)  do  not  include  the  noise  or 
measurement  errors.  However,  the  errors  of  the  observed  points  can 
be  included  in  the  predictions  by  replacing  in  equations  (9)  and 
(12)  with 

fik  ° cik  * Dik’ 


Dik  ‘ M(ni"k1 


is  the  autocovariance  of  the  noise. 


The  error  covariance  is  the  correlation  of  the  predic- 


tion errors  ep  and  at  points  P and  Q, 


L 


1 

I 

B 


2 


a 


> 


► 


. 


°PQ  = M{epeQ}  = CpQ 


n 


kl}  CPk  Cik_1  CQi 


n’  n' 


' Ji  c«*  V CpJ 


n n n'  n' 


+ y y y y cD.  c..-1  c. . c.  -1  cn  (is) 

1=1  k=i  j=i  t=i  pk  ,k  ’j  " 


where  the  observed  points  (j,  a)  associated  with  point  Q do  not 
necessarily  coincide  with  the  (i,  k)  points.  The  relationship  in 
equation  (11)  is  a special  case  with  P=Q. 

B.  Profile  Prediction  Method 

Profile  data  may  also  be  regarded  as  sample  points  of 
a continuous  function.  Then  the  methods  of  spectral  analysis  are 
available  for  the  prediction  of  the  gridded  data  base.  The  geoid 
height  autocovariance  function  of  Tscheming  and  Rapp  [4]  is  used  in 
the  prediction  of  the  mean  geoid  heights  and  the  prediction  error 
covariance  function.  The  Fourier  transform  of  the  geoid  height 
covariance  function  is  the  power  spectrum  of  the  geoid  height  field 
which  is  assumed  to  be  isotropic.  The  transform  changes  the  data  on 
a flat  surface  from  the  time  or  space  domain  to  the  frequency  domain. 
The  Wiener  theory  of  prediction  is  used  to  estimate  the  mean  geoid 
heights  and  the  prediction  error  spectrum  in  the  frequency  domain. 
Then  the  inverse  Fourier  transform  of  the  error  spectrum  is  the 
prediction  error  autocovariance.  The  profile  prediction  method  is 
discussed  in  the  appendix. 


• • 
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IV.  ESTIMATION  OF  GRAVITY  ANOMALIES 

The  gravity  anomaly  estimates  presented  in  this  report  are 
computed  from  simulated  geoid  heights  and  are  compared  with  gravity 
anomalies  known  to  be  consistent  with  the  input  geoid  heights.  This 
section  introduces  Stokes'  integral  and  a method  of  localization. 

The  localization  is  supported  by  a discussion  of  the  contributions  of 
distant  zones  and  high  degree  harmonic  coefficients  to  the  computed 
geoid  height.  Effects  of  errors  in  the  low  degree  harmonic  coeffi- 
cients of  the  earth  gravitational  model  (EGM)  are  also  considered. 

A.  Stokes'  Integral 

The  geoid  height  N and  the  gravity  anomaly  Ag  are 
related  to  the  disturbing  potential 
T = W-U, 

where  W is  the  total  gravity  potential  and  U is  the  normal  gravity 
potential  (the  potential  field  of  the  reference  ellipsoid).  The 
disturbing  potential  is  representable  as  a harmonic  series 

00  n 

TU’X)  = ^ n=2  m=0  CCnm  “S  mX  + S"”  Sln  "U)Pn>»<Sln*):i-  (,6) 

(n,m)  1 (2,0),  (4,0). 

The  gravity  anomaly  is  defined  in  terms  of  the  disturbing  potential 


£L  . 21 

ar  ~ r 


and  the  geoid  height  is  defined  by  Bruns'  equation 


N = — . 
Y 
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Stokes'  Integral  makes  it  possible  to  compute  N from  Ag  without 
working  through  the  disturbing  potential.  Stokes'  Integral  is 


N(*,x)  = / / S (ip) Ag ( <t> ' ,x'  ) cos*1  d*'  dx*  (17) 

over  an  approximating  sphere  of  mean  radius  R and  mean  gravity  y. 
Stokes'  function  for  the  angular  distance  * between  (<|>,x)  and 
(<t>V)  is 

S(*)  = — ““t 6 sin(^)  + 1-5  cos* 

slntf)  2 

- 3 cos*  tn[sin(|r)  + sin2(|-)]  (18) 

for  points  on  the  unit  sphere  [5]. 

Stokes'  integral  is  ordinarily  used  to  compute  geoid 
heights  by  worldwide  integration  of  the  gravity  anomaly  field.  For 
this  purpose  the  integral  is  replaced  by  the  discrete  form 

N(*,x)  = l l S(*)Ag(*',x')  cos*'  a*'  ax'  (19) 

4 Y *'  x' 

with  some  provision  for  dealing  with  the  singularity  of  S(*)  at  *=0, 
and  Ag  Is  the  mean  gravity  anomaly  over  an  element  with  center  ( * ' ,X ' ) . 

Stokes'  integral  may  also  be  regarded  as  an  integral 
equation  for  the  inverse  problem  of  estimating  anomalies  from  geoid 
heights.  Then  the  discrete  form  (19)  becomes  a worldwide  linear 
algebraic  system 
Ax  = b, 

where  A is  the  coefficient  matrix  of  values  of  Stokes'  function,  b 
is  the  data  vector  of  geoid  height  values,  and  x is  the  system 


parameter  vector  of  gravity  anomalies.  Since  reliable  geoid  height 
values  are  not  available  for  much  of  the  world  and  the  worldwide 
algebraic  system  is  much  too  large  to  estimate  the  vector  x over  a 
one  degree  grid,  a method  of  localization  of  the  computations  is 
needed.  That  is,  there  is  a need  to  restrict  the  data  and  computa- 
tions for  Stokes'  integral  equation  to  a relatively  small  local 
region.  The  next  section  proposes  a method  of  localization  which 
consists  of  estimating  the  global  (long  wave  length)  harmonics  with 
the  earth  gravitational  model  of  a world  geodetic  system,  and  then 
estimating  the  high  degree  (short  wave  length)  components  over  a 
small  local  region. 

B.  Localization  of  Stokes'  Integral  Equation 

Each  of  the  functions  N,  S(^),  and  Ag  involved  in 
Stokes'  integral  equation  is  harmonic  and  expandable  in  a series  of 
Legendre  associated  functions.  The  harmonic  series  are 

oo  p 

N(<m)  = 1 l C(Cnm  cos  mx  + Snm  sin  mx)Pnrn(sin*)] , (20) 

SM  = l vj-  pn(cos*)’  (21) 

n=2  n_l  n 

Pp (costji)  = Pn (cos«)  Pn(cos<j>') 


oo  n 


AgU.x)  = — I l (n-1  )[C  cos  mx 
R2  n=2  m=0  nm 


+ Snm  sin  mx]Pnm(sin*),  (23) 


where  (n,m)  f (2,0),  (4,0). 


Since  the  Legendre  functions  are  orthogonal,  they  are  linearly  inde- 
pendent so  that  N,  S(i ii) , and  Ag  may  be  decomposed  into  high  and  low 
degree  components.  Figure  5 shows  the  effect  of  removing  lower  har- 
monics from  Stokes'  function,  with  the  full  function  and  three 
residual  forms  being  shown.  The  high  degree  (short  wave  length) 
components  have  short  distance  effects  and  the  low  degree  components 
have  long  range  effects. 

Linearity  of  the  integral  allows  a corresponding  decom- 
position of  the  surface  of  the  earth  according  to  distance  from  the 
point  of  estimation.  Then  the  worldwide  integral  is  the  sum  of  the 
integrals  over  a near  region  R-|  and  a distant  region  R^  such  that 
the  tv/o  regions  together  cover  the  surface  of  the  earth.  In  this 
report  it  is  assumed  that  the  low  degree  component  of  Stokes'  func- 
tion consists  of  harmonic  degrees  2-12  and  the  high  degree  component 
consists  of  degrees  13-®.  It  is  further  assumed  that  the  region  R1 
has  a radius  not  less  than  5°  for  1°  x 1°  estimation. 

By  means  of  the  decompositions  mentioned  above,  Stokes' 
integral  equation  may  be  written 


N(<M)  = K / R / S2_12(i^)Ag2_12(4.' , a'  ) cos <t> ' d $'  dx' 


18 


Q 

B 


I 


» 


+ K/  / S13_oo(^)Ag13_oo(<t>',X,)cost'  d<j>'  dx* 

R1 

+ K / / s2_-j 2 ((p)Ag2_-1 2 (<#) 1 ,A‘  )cos0‘  <4'  dA ' 

r2 

+ K / / 3_0o( ^ ) Ag-j 3_00 ( ‘ ' ) cos (+> ' cV  dx,  (24) 

r2 

where  K = . 

4iry 

The  first  and  third  integrals  on  the  right  side  of 

equation  (24)  are  global  components  and  may  be  approximated  with 

harmonic  series  in  C , S of  the  earth  gravitational  model.  With 

nm  nm 

suitable  choice  of  R1  and  R2>  the  fourth  integral  may  be  considered 
negligible  due  to  the  high  harmonic  degree  of  the  integral  and  the 
distance  of  the  region  from  the  computation  point.  Thus,  Stokes' 
integral  equation  reduces  to 

Ni3_.  “ 477  / / S13_>)^13_jV,X')coS*'  d*'  dx\  (25) 

R1 

by  removing  the  low  degree  components  of  the  geoid  heights;  that  is, 

'^13_oc  = ^2-«>  " ^2-12  (26) 

The  approximation  (25)  allows  computation  of  the  high  degree  estimate 
Agp_oo  in  a small  local  region.  With  the  high  degree  estimate  Ag.^^, 
the  total  anomaly  estimate  Ag^  is 

“9 2 -«  = ig2-12  + Agl3-oo* 

Methods  of  estimation  of  Ag-j^  are  discussed  in  more  detail  in 
Section  IV. E. 
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C.  Contributions  of  Distant  Zones  and  High  Harmonic  Degrees 
To  Computed  Geoid  Height 

The  expected  error  variance  o2(N)  of  the  computed  geoid 
height  due  to  neglect  of  gravity  anomaly  data  beyond  a region  of 
radius  and  harmonic  degree  greater  than  K is,  from  Heiskanen  and 
Moritz  [5], 

”2(N)  ’ $r  „X,  Q2n  Cn’  <28> 

where 

7T 

Q„  * / S(g,)Pn(cos»)sin»  d*,  (29) 

♦o 

Cn  is  the  degree  variance  of  degree  n, 

S ( <p)  is  Stokes'  function  of  angular  distance  ip. 

Computed  results  for  the  degree  variances  of  the  test  model  are  shown 
graphically  in  Figure  6.  The  graph  shows  the  expected  error  in  geoid 
heights  due  to  neglecting  gravity  anomaly  data  outside  a circle  of 
radius  <1^  and  harmonics  above  n = 8,  12,  and  18.  For  high  degree 
components  above  degree  12,  the  graph  shows  that  the  expected  error 
for  i^g  = 5°  is  approximately  1.5  meters. 

The  inverse  problem  of  estimating  the  expected  error 
variance  a2(Ag)  of  gravity  anomalies  due  to  neglect  of  high  degree 
harmonics  and  distant  regions  is  much  more  difficult.  The  error 
variance  for  geoid  heights  suggests  that  high  degree  gravity  anomaly 
components  of  degree  greater  than  12  may  be  effectively  estimated 
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Figure  6.  Expected  Error  in  N Due 


by  neglecting  geoid  height  data  at  distances  greater  than  5°.  It 


L 


is  assumed  that  the  fourth  integral,  over  Rg,  of  equation  (24)  is 
negligible  since  the  high  degree  Stokes'  function  is  effectively 
damped  out  for  the  distances  included  in  this  integral  (Figure  5). 
The  present  series  of  tests  is  a more  specific  justification  for 
this  assumption  than  the  statistical  method  discussed  above. 

D.  Effects  of  Errors  in  the  Low  Degree  Component 

Errors  in  the  low  degree  harmonic  coefficients  are  not 
related  to  the  satellite  altimetry  data,  but  are  included  in  the 
computational  process  by  means  of  the  supplementary  data.  They  have 
a two-fold  effect  on  the  estimated  gravity  anomalies,  which  can  be 
seen  by  analyzing  the  computational  process.  They  affect  the  high 
degree  gravity  anomaly  estimates  by  being  included  in  the  geoid 
height  data  used  in  the  recovery  process,  and  also  they  have  an 
effect  on  the  low  degree  component  since  it  is  computed  from  the 
coefficients.  An  estimate  of  the  magnitude  of  this  effect  was 
determined  by  repeating  the  recovery  process  with  perturbed  coeffi- 
cients. The  results  of  this  analysis  show  that  the  effect  on  the 
geoid  height  values  is  ±0.6  meters,  which  causes  an  error  of  ±0.2 
mgal  in  the  high  degree  component  of  the  gravity  anomalies.  The 
effect  of  these  errors  on  the  low  degree  components  is  about  ±0.7 
mgal,  and  therefore  the  combined  effect  on  the  gravity  anomalies 
is  ±0.7  mgal.  These  errors  in  the  harmonic  coefficients  will  be 
reduced  by  updating  the  earth  gravitational  model  as  geoid  height 
data  from  satellite  altimetry  becomes  available. 
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E.  Method  of  Estimating  Anomalies 

The  discrete  form  of  N-j 3_00(R-| ) » equation  (25),  may  be 
solved  as  a linear  algebraic  system 

Ax  = b,  (30) 

where  A is  the  coefficient  matrix  with  one  row  and  one  column  for  each 
square  in  region  Rj.  For  a least  squares  estimator,  R1  may  be  decom- 
posed into  an  inner  and  outer  region  with 

Ax  + v = b,  (31) 

where  v is  the  residual  vector  and  A has  one  row  for  each  square  in 
the  outer  region  and  one  column  for  each  square  in  the  inner  region. 

With  the  least  squares  estimator,  anomalies  are  estimated 
only  In  the  inner  region  and  there  is  a loss  of  information  in  the 
outer  region.  Both  methods  of  estimation  lose  accuracy  near  the 
boundary.  These  errors  of  localization  affect  the  parameter  vector 
x and  the  coefficient  matrix  A so  that  the  errors  carry  over  to  all 
elements  of  the  covariance  matrix.  Thus,  even  though  only  the  esti- 
mate of  the  central  square  is  accepted  for  each  run,  the  corresponding 
element  of  the  covariance  matrix  is  dominated  by  localization  error 
and  the  estimation  error  must  be  evaluated  by  some  other  method. 

Most  of  the  test  runs  presented  in  this  report  were  made 
with  an  unweighted  least  squares  estimator.  Then  the  estimator  x of 
x for  equation  (31 ) is 

x = (ATA)~ 1 ATb.  (32) 

The  outer  region  was  11°  x 11°  and  the  inner  region  was  5°  x 5°.  The 
singularity  of  S ( ip)  was  avoided  by  arbitrarily  offsetting  by  eight 
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minutes  from  ^ = 0.  The  program  was  modified  for  iterative  improve- 
ment of  two  estimated  anomaly  fields.  Estimation  errors  were  com- 
puted by  subtracting  anomalies  known  to  be  consistent  with  the  input 
geoid  heights  from  the  estimated  anomalies. 

In  these  computations,  the  eight  minute  offset  of  S (0) 
was  used  to  avoid  the  singularity  of  Stokes'  function.  However, 
additional  tests  suggest  that  a larger  offset  may  have  been  better 
for  the  1°  x 1°  estimates.  An  optimum  offset  of  18'  was  selected 
after  analyzing  offset  values  out  to  about  21'  from  the  computation 
point.  Moritz  [5]  discusses  a method  in  which  the  size  of  the  off- 
set depends  on  both  the  size  of  the  square  and  the  probable  error 
of  the  geoid  field.  Further  study  of  this  approach  is  being  carried 
on  at  this  Center. 

It  is  also  possible  to  introduce  a local  correction  to 
the  harmonic  coefficient  Cqq  or  the  equatorial  radius  a.  This  might 
be  done  by  subtracting  the  mean  geoid  height  over  the  integration 
region  from  the  geoid  height  of  each  one  degree  square  or  by  inserting 
a column  of  ones  into  the  design  matrix  A of  the  least  squares  estima- 
tor to  estimate  the  local  bias.  The  use  of  weight  matrices  has  not 
been  completely  tested.  The  best  weight  matrix  for  this  purpose 
appears  to  be  the  inverse  of  the  geoid  height  covariance  matrix.  A 
least  squares  estimator  weighted  in  this  way  becomes  a minimum  vari- 
ance estimator.  Another  procedure  that  will  be  investigated  is  the  use 
of  initial  estimates  of  the  gravity  anomalies  based  on  existing  gravity 
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observations.  The  computer  program  could  be  modified  for  any  of 

these  methods.  A computer  program  for  least  squares  collocation  [6]  • • 

is  being  prepared  to  continue  these  lines  of  investigation. 
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V. 


TEST  DATA  SIMULATION 


Tests  were  conducted  to  demonstrate  that  gravity  anomalies 
can  be  predicted  from  local  geoid  heights  with  sufficient  accuracy 
to  make  effective  use  of  GEOS- III  altimetric  geoid  height  data. 

Testing  of  the  computational  method  described  earlier  requires  con- 
sistent sets  of  simulated  high  frequency  geoid  heights  and  gravity 
anomalies.  In  some  of  the  tests,  the  geoid  heights  were  modified 
by  the  introduction  of  random  errors  or  by  smoothing  to  simulate  mean 
values. 

Three  methods  of  developing  consistent  sets  of  high  fre- 
quency geoid  heights  and  gravity  anomalies  are  considered  feasible: 
the  gravity  anomaly  method,  the  point  mass  method,  and  the  harmonic 
method.  In  the  first  method,  a gravity  anomaly  field  is  assumed  and 
the  geoid  heights  are  computed  with  Stokes'  integral  and  worldwide 
integration.  In  the  second  method,  a point  mass  set  replaces  the 
mass  of  the  earth  in  the  computations  and  both  the  gravity  anomaly 
field  and  the  geoid  height  field  are  computed  from  the  disturbing 
potential  field  of  the  point  mass  set.  In  the  third  method,  a set 
of  harmonic  coefficients  to  high  degree  and  order  is  constructed 
and  both  the  gravity  anomaly  and  geoid  height  fields  are  computed 
by  harmonic  series.  The  last  method  appears  to  be  the  simplest  and 
most  consistent  since  no  worldwide  integration  is  involved.  Moreover, 
the  exact  harmonic  structure  of  the  data  is  known  so  that  separation 
of  the  high  and  low  degree  components  of  geoid  heights  or  gravity 
anomalies  is  simple. 
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In  order  to  produce  a gridded  data  base  whose  exact  harmonic 


content  is  known  for  testing  the  estimation  of  gravity  anomalies,  it 
was  decided  to  generate  the  data  base  from  a hypothetical  harmonic 
model  consisting  of  spherical  harmonic  coefficients  to  degree  and 
order  18  from  an  existing  earth  gravitational  model  and  fictitious 
coefficients  beginning  from  degree  19  through  degree  and  order  141. 
The  fictitious  coefficients  were  generated  by  Kaula's  approximate 
formula  [7] 

C = — — mgal2  (33 

n n/i  i 3'  3 

n 1+2rT 

and  the  harmonic  form 


C = M{ Ag2  } = V (A  2 + B 2)  ( 

n n ' nm  nm 

so  that  all  A and  Ef  of  degree  n have  the  same  magnitude  and  random 

signs,  and  the  Cn  in  equations  (33)  and  (34)  agree.  The  l\m>  Bnm 

were  transformed  to  C , S which  were  used  to  compute  compatible 

nm  nm  r 

sets  of  geoid  heights  and  gravity  anomalies  at  center  points  of 
1°  x 1°  elements  of  a 50°  x 50°  block  with 


141  n 


N<*.»>  - l l (c  cos  « ♦ Snm  sin  «)P  (sin*),  (35) 

1 n -L  m=0 


&g(<M)  =—  l (n-1 ) l (C  cos  mx 
R2  n=2  m=0  nm 


+ Snm  sin  mX)Pnm(sin<j>) . (36) 


A 10°  x 10°  block  located  in  the  center  of  the  data  was  selected  for 
testing  and  the  geoid  heights  were  input  to  the  estimation  program 
which  subtracted  the  low  degree  component  and  then  estimated  high 
degree  anomalies  from  the  high  degree  geoid  heights.  The  low  degree 
component  Ag(2-12)  of  the  anomaly  field  was  added  to  the  estimated 
Ag ( 1 3-1 41 ) and  the  results  compared  with  the  simulated  Ag (2-1 41 ) . 
Although  the  simulated  gravity  anomalies  are  not  intended  to  repre- 
sent the  physical  values,  they  are  consistent  with  the  simulated 
geoid  heights.  Hence,  the  differences  between  the  estimated  gravity 
anomalies  and  the  simulated  values  are  the  errors  of  estimation. 
Figure  7 is  a flow  chart  of  the  test  procedure  for  the  hypothetical 


VI.  PROGRAM  INPUT  AND  TEST  RESULTS 

Results  from  the  estimation  of  gravity  anomalies  of  one 
hundred  1°  x 1°  squares  in  a 10°  x 10°  block  from  simulated  input 
data  and  the  associated  estimation  errors  are  displayed  in  tabular 
and  graphic  form.  The  graphics  include  contour  plots  of  the  data 
and  amplitude  spectra.  Mean  gravity  anomaly  values  were  simulated 
in  these  computations  by  a smoothing  operation.  The  effect  of 
B measurement  error  in  the  geoid  height  values  was  approximated  by 

adding  random  noise  to  the  data. 

Contour  plots  of  the  geoid  heights  and  gravity  anomalies 
, at  the  centers  of  1°  x 1°  squares  for  the  full  spherical  harmonic 

set  (n,m  = 2-141)  are  shown  in  Figures  8 and  9,  respectively.  These 
fields  are  considered  errorless.  Point  gravity  anomalies  from  the 
g high  degree  component  (n,m  = 13-141)  of  the  anomaly  field  for  each 

1°  x 1°  element  are  shown  in  Table  1.  A set  of  1°  x 1°  mean  gravity 
anomalies  was  simulated  by  applying  two-dimensional  smoothing  to 
the  high  degree  anomaly  field.  The  smoothed  value  for  each  element 
is  one-half  ':ts  raw  value  plus  one-eighth  of  the  sum  of  the  values 
of  the  four  adjacent  squares.  The  smoothed  field  is  shown  in 
Table  2. 

The  gravity  anomalies  estimated  from  gridded  geoid  heights 
for  1°  x 1°  elements  are  shown  in  Table  3.  To  determine  the  estima- 
, tion  errors  for  each  1°  x 1°  element,  differences  were  computed 

between  the  gravity  anomaly  estimates  (Table  3)  and  both  the  reference 
point  values  of  Table  1 and  the  reference  mean  values  of  Table  2. 
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Table  1 


1°  x 1°  POINT  GRAVITY  ANOMALIES  FROM  THE 
HYPOTHETICAL  HIGH  DEGREE  FIELD  n,m  = 13+141 
(Ag,  mgal) 


-21.5 

12.7 

38.1 

33.7 

26.6 

36.6 

46.4 

39.9 

25.4 

13.6 

- 7.1 

- 3.2 

20.7 

37.8 

32.1 

24.5 

29.4 

26.7 

n 

-15.9 

4.0 

- 8.0 

3.6 

29.5 

33.6 

25.4 

25.6 

17.7 

- 6.3 

-21.1 

- 3.1 

- 8.3 

3.5 

28.0 

27.5 

13.4 

11.0 

H 

- 3.3 

-10.9 

- 0.8 

-22.3 

-11.6 

20.0 

22.9 

9.6 

12.1 

21 .4 

24.1 

15.6 

13.9 

-22.4 

-22.2 

n 

13.7 

■ 

6.8 

15.6 

18.5 

11.8 

n 

-14.6 

-12.3 

H 

H 

- 6.0 

-11.2 

-12.3 

-14.0 

-18.5 

-27.0 

-32.7 

-16.0 

0.6 

8.0 

12.0 

12.1 

10.2 

9.6 

2.0 

-10.4 

-14.0 

1.7 

10.5 

15.5 

25.7 

28.4 

22.7 

19.5 

12.6 

23.5 

23.9 

22.1 

13.2 

..  ---° 

13.7 

13.3 

6.9 

4.8 

3.9 

MEAN  VALUE 
RMS 


8.7  mgal 
19.0  mgal 


Table  2 
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1°  x 1°  SMOOTHED  MEAN  GRAVITY  ANOMALIES  FROM 
THE  HYPOTHETICAL  HIGH  DEGREE  FIELD  n,m  = 13*141 
(Ag,  mgal) 


■15.3  11.0  31.1  31.1  28.1  30.8  41.6  34.0  19.9  9.1 


- 6.6  0.7  19.9  33.4  31.4  27.7  30.1  24.8  6.0  -10.8 


BH 


- 1.1  -16.5  - 8.4  15.8  20.3  11.4  12.2  18.1  18.6  10.3 


— 
31 .1 

— 

31.1 

19.9 

33.4 

n 

27.6 

3.2 

24.1 

- 8.4 

15.8 

-16.0 

5.2 

-12.3 

1 .6 

-13.3 

1.0 

1.2 

9.1 

14.9 

10.6 

8.6  -16.8  -16.0  5.2  11.6  5.2  6.0  12.1  13.9 


■ 


30.8 

41.6 

27.7 

30.1 

24.8 

23.2 

15.9 

12.8 

11.4 

12.2 

5.2 

6.0 

- 1.8 

- 5.5 

11.0 

11.0 

21 .6 

23.4 

13.4 

13.3 

19.9 

9.1 

6.0 

-10.8 

n 

-12.7 

0.2 

- 8.3 

18.6 

10.3 

■ 


■ 


.3  -12.2 


18.8  -25.3  -13.3  1.0  8.0  11.0  11.0  9.1  7.0  - 0.4 


- 6.5  - 9.2  1.2  9.1  14.4  21.6  23.4  19.5  16.0  8.7 


1H 


MEAN  VALUE  = 8.3  mgal 

RMS  = 16.3  mgal 
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For  this  evaluation,  Tables  1 and  2 are  considered  errorless.  The 
differences  are  shown  in  Tables  4 and  5,  respectively.  The  root 
mean  square  (RMS)  error  of  estimates  from  point  data  is  6.1  mgal  and 
the  RMS  error  of  estimates  from  mean  data  is  2.9  mgal.  Thus,  1°  x 1° 
estimates  are  better  in  comparison  with  mean  data  than  with  point 
data. 

Input  and  output  data  fields  were  analyzed  by  spectral 
methods  for  their  harmonic  frequency  content.  Contour  plots  of  the 
amplitude  spectra  are  shown  in  Figures  10  to  14.  The  sampling  inter- 
val on  the  contour  plots  is  a one  degree  increment  of  latitude  or 
longitude  between  data  points.  The  coordinate  axes  are  frequencies 
of  the  Fourier  representation  of  gravity  anomalies  or  estimation 
errors.  The  wave  lengths  in  degrees  are  reciprocals  of  the  frequen- 
cies. The  axes  are  in  the  direction  of  increasing  frequency  and 
decreasing  wave  length. 

The  amplitude  spectra  of  gravity  anomalies  and  errors  were 
computed  as  though  the  data  points  were  on  a plane.  Contours  of 
amplitude  spectra  for  isotropic  fields  over  a plane  are  arcs  of 
concentric  circles.  Curvature  of  the  surface  and  convergence  of 
meridians  distort  the  contours  of  an  isotropic  field  to  elliptical 
arcs.  If  the  field  is  not  isotropic,  then  the  contours  are  dis- 
torted accordingly. 

The  contour  plots  of  amplitude  spectra  in  Figures  10,  11 
and  12  show  that  the  fields  are  approximately  isotropic  in  the  long 
wave  lengths.  There  is  distortion  in  the  contours  of  the  high 
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Table  4 

ERROR  OF  ESTIMATION  (&~g  - Ag) 
FOR  1°  x 1°  POINT  VALUES 
(mgal) 


5°  x 5°  INTEGRATION 
MEAN  VALUE  = -0.2  mgal 

RMS  = 6.1  mgal 


Table  5 


ERROR  OF  ESTIMATION  (A~g  - Ag) 
FOR  1°  x 1°  MEAN  VALUES 
(mgal) 


• • 


1.1  -1.2  - 3.2  - 3.1  - 2.7  0.8  - 5.5 


0.7  0.2  - 2.1  - 3.6  - 3.3  - 3. 


n 


IB 


1.2  4.1  2.4  - 2.4  - 2.8  - 0.3 


6.5  9.2  6 


2.5  0.8  0.3 


2.0  - 1.3  - 0.2 


5°  x 5°  INTEGRATION 
MEAN  VALUE  - 0.2  mgal 

RMS  = 2.9  mgal 


B 

- 2.2 

- 1.0 

- 3.6 

- 0.7 

1.6 

- 1.7 

1.3 

- 0.5 

- 0.6 

0.3 

1.4 

- 2.6 

- 3.8 

- 2.9 

- 1.3 

- 3.3 

- 2.9 

B 

2.3 

2.6 

0 

- 0.1 

0.9 

- 3.1 

- 2.6 

- 0.9 

1.5 

1.4 

1.3 

• • 
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Fiquro  10.  Amolitude  Spectrum  of  1°  x 1’  Point  Gravity  'Inc,  alios 


ro 


Fiqure  11.  Amplitude  Spectrum  of  1°  x 1°  Mean  Gravity  Anomalies  (Aq) 
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of  latitude 
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Figure  12.  Amplitude  Spectrum  of  1°  x 1°  Estimated  Gravity  Anomalies  (A~g). 
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Figure  13.  Amplitude  Spectrum  of  1°  x 1°  Estimation  Errors  (Ag-Ag), 
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Figure  14.  Amplitude  Spectrum  of  1°  x 1°  Estimation  Errors  (A~g-Aq). 


I 


frequency  components.  Comparison  of  spectra  for  point  gravity 
anomalies  (Ag ) and  mean  (smoothed)  gravity  anomalies  ( Ag”)  shows  that 
there  is  a loss  of  high  frequency  in  smoothing.  Figures  11  and  12 
are  almost  identical.  This  indicates  that  the  estimated  gravity 
anomaly  field  and  the  smoothed  gravity  anomaly  field  have  very  nearly 
the  same  Fourier  harmonic  content.  The  closer  similarity  of  ampli- 
tude spectra  for  smoothed  values  and  estimated  values  is  corroborated 
by  the  spectra  for  estimation  errors  in  Figures  13  (estimated  minus 
point  values)  and  14  (estimated  minus  mean). 

The  preceding  results  were  obtained  with  an  errorless  geoid 
height  field.  Since  measurement  necessarily  involves  error,  the 
effects  of  random  error  and  bias  in  the  geoid  field  were  investigated 
and  the  results  shown  in  Tables  6 and  7.  A random  error  field  with  a 
distribution  of  zero  mean  (y=0)  and  standard  deviation  (o)  of  one 
meter  was  added  to  the  geoid  field.  The  error  of  estimation  for  each 
1°  is  the  difference  between  the  resulting  predicted  anomalies  (Ag) 
and  the  simulated  mean  anomalies  (Table  6).  Inclusion  of  the  p=0, 
a=lm  errors  in  the  geoid  field  increases  the  RMS  anomaly  error  from 
2.9  mgal  to  9.0  mgal  with  no  significant  change  in  the  mean  error. 

A one  meter  bias  was  then  added  to  the  geoid  field.  The  error  of 

A 

estimation  for  this  set  of  anomalies  (Ag)  is  shown  in  Table  7.  Addi- 
tion of  the  one  meter  bias  increased  the  mean  error  from  0.2  mgal  to 
1.5  mgal  and  the  RMS  from  2.9  mgal  to  3.5  mgal. 

Tests  for  the  effects  of  random  error  in  the  geoid  field 
were  extended  to  include  random  errors  with  o = 0.3,  0.6,  1.2,  and 


1 
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Table  6 

ERROR  OF  ESTIMATION  (Ag-Ag)  FOR  l°xl°  MEAN  VALUES 
FROM  A GEO ID  WITH  A RANDOM  ERROR  OF  ONE  METER 


15.9 

- 9.3 

-16.9 

-26.3 

17.7 

- 9.3 

- 3.1 

-11.7 

-12.2 

- 6.3 

m 

6.9 

2.6 

-11.0 

-21.7 

- 3.8 

- 2.8 

0 

n 

- 8.3 

- 0.7 

- 0.6 

5.6 

16.5 

-10.3 

■ 

0.5 

14.8 

n 

2.9 

m 

3.9 

- 1.2 

- 3.9 

- 1.1 

-12.9 

-11.7 

15.1 

- 1.3 

H 

16.3 

- 1.8 

- 0.5 

6.1 

13.9 

-11.7 

8.4 

- 6.6 

- 2.8 

4.0 

11.4 

2.9 

19 

1 .9 

-14.2 

4.1 

2.2 

- 2.5 

9.3 

- 5.6 

1.3 

- 0.3 

■ 

- 5.0 

3.3 

3.8 

21.2 

0.9 

■ 

22.7 

14.1 

n 

6.2 

H 

0.6 

m 

- 8.7 

10.0 

0.8 

H 

5.3 

- 5.3 

-12.6 

-10.9 

3.6 

- 8.0 

n 

5.9 

- 5.2 

- 1.9 

- 1.2 

18.1 

5.3 

- 3.2 

- 6.2 

- 1.6 

■ 

5°  x 5°  INTEGRATION 
MEAN  VALUE  = 0.1  mgal 

RMS  = 9.0  mgal 


MEAN  VALUE  = -1.5mgal 


RMS 


3.5  mgal 
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3.0  meters,  all  with  zero  mean.  The  error  of  estimation  increased 
from  2.9  mgal  to  26.2  mgal  as  a increased  from  0 to  3 meters.  Results 
of  this  analysis  are  shown  in  Figure  15.  Also  shown  are  results  for 
2.5°  and  5°  elements  accomplished  by  averaging  the  1°  results,  and, 
as  expected,  the  error  is  reduced  as  the  element  size  is  increased. 

For  the  purpose  of  iterative  improvement,  new  geoid  heights 
were  computed  from  the  estimated  anomalies  by  Stokes'  integral  and 
a &N  field  was  obtained  by  subtracting  the  true  geoid  heights  from 
the  new  geoid  heights.  The  AN  values  were  input  to  the  estimation 
program  and  corrections  to  the  original  anomaly  estimates  were  com- 
puted with  an  algebraic  solution  over  a 5°  x 5°  integration  region. 

The  corrected  estimates  are  compared  with  point  and  mean  anomalies 
in  Table  8. 

Table  8 

ITERATIVE  IMPROVEMENT 
(36  1°  x 1°  SQUARES) 


RMS 

MEAN 

ANOMALIES 

POINT 

ANOMALIES 

First  Estimate 
Iteration 

2.20  mgal 
1 .45  mgal 

4.42  mgal 
2.23  mgal 

(meters) 
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VII.  CONCLUSIONS 

It  is  possible  to  localize  the  estimation  of  gravity 
anomalies  from  a 1°  x 1°  gridded  data  base  of  geoid  height  values. 

The  method  of  localization  consists  of: 

Approximating  the  low  degree  (long  wave  length)  component 
of  gravity  anomalies  and  geoid  heights  from  the  harmonic  coefficients 
of  an  earth  gravitational  model. 

Estimating  the  high  degree  (short  wave  length)  components 
of  gravity  anomalies  from  high  degree  components  of  geoid  heights 
over  a local  region  with  a linear  estimator. 

Mean  gravity  anomalies  for  1°  x 1°  squares  can  be  effec- 
tively estimated  from  localized  geoid  heights  which  contain  random 
errors  whose  standard  deviation  is  one  meter  or  less.  Random  errors 
with  a standard  deviation  of  three  meters  increase  the  RMS  error  of 
estimation  to  26  mgal.  Mean  anomalies  for  5°  x 5°  squares  can  be 
locally  estimated  from  1°  x 1°  geoid  heights  whose  random  errors 
have  a standard  deviation  as  large  as  three  meters.  A bias  of  one 
meter  in  the  local  geoid  height  field,  if  not  removed,  has  a signi- 
ficant effect  on  the  mean  error  of  the  gravity  anomaly  estimates, 
although  such  a bias  does  not  change  the  RMS  error  of  estimation 
appreciably. 

Errors  in  the  low  degree  harmonic  coefficients  of  the 
EGM  affect  both  the  high  and  low  degree  components  of  the  anomaly 
estimates.  These  errors  can  be  reduced  by  updating  the  EGM  as 
altimetry  data  accumulates.  Methods  of  avoiding  the  singularity  of 
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Stokes’  kernel  function,  the  use  of  smaller  element  sizes  and  methods 

of  minimizing  the  effects  of  bias  in  the  data  are  indicated  to  be  * 

areas  for  further  study. 

The  unweighted  least  squares  estimator  is  not  necessarily 
the  best  choice  of  estimation  model  since  it  does  not  use  all  avail-  ' * 

able  data.  The  model  could  be  modified  to  use  weight  matrices 
derived  from  covariance  matrices  and  to  compute  corrections  for 

initial  estimates  of  gravity  anomalies.  Moreover,  the  selected  * ® 

estimator  apparently  does  not  achieve  the  optimum  resolution  of 
the  data  into  signal  (true  values)  and  noise  (random  errors). 

Additional  study  is  also  indicated  in  this  area. 

»•  • 
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APPENDIX 
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PREDICTION  OF  GEOID  HEIGHTS  FROM  PROFILE  DATA 


Al.  INTRODUCTION 


Profile  data  for  GEOS- I I I altimetry  consists  of  measured 
geoid  heights  at  closely  spaced  intervals  along  approximately  parallel 
tracks.  This  appendix  outlines  a method  for  estimating  the  predic- 
tion error  covariance  function  for  Wiener  prediction  of  a geoid  height 
gridded  data  base  from  profile  data.  The  error  covariance  function 
is  computed  in  the  frequency  domain  for  simplicity  since  integrations 
in  the  time  domain  become  multiplications  in  the  frequency  domain. 

The  accuracy  of  the  predicted  1°  x 1°  geoid  height  values  is  given  by 
the  error  integral,  which  is  the  integral  of  the  error  covariance 
function  over  the  1°  x 1°  surface  elements.  Also  included  in  the 
discussion  are  geoid  height  autocovariance,  power  spectra,  geoid 
height  prediction,  the  spectral  coefficients  of  the  Wiener  equation 
and  their  application  in  the  error  analysis  of  predicted  geoid 
heights. 

A2.  THE  GEOID  HEIGHT  AUTOCOVARIANCE  FUNCTION 


The  prediction  method  is  based  on  the  assumption  that  the 
geoid  heiqht  field  can  be  treated  as  an  isotropic  stationary  stochastic 
process.  The  effect  of  isotropy  for  a covariance  function  between 
two  profiles  with  an  origin  at  point  P is  shown  in  Figure  Al . The 
value  of  the  geoid  height  autocovariance  function  relating  points  P 
and  Q depends  only  on  the  angular  distance  between  P and  Q and  not 
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on  the  location  of  the  points  or  the  direction  of  one  from  the  other. 
The  line  segment  PQ ' is  normal  to  both  profiles. 


Figure  A1 . Geoid  Height  Autocovariance  Function 


In  these  computations,  a form  of  the  geoid  height  autocovariance 
function  of  Tscherning  and  Rapp  [4],  with  the  low  degree  and  the 
high  degree  parts  removed,  is  used: 


CM 


1080  C 
n=13  ( n- 1 ) 2 


Fn+1  Pn(cos*), 


(Al) 


where 


R = Mean  radius  of  earth, 

y = Theoretical  gravity  of  the  reference  ellipsoid, 
Pn(cosijj)  = Legendre  polynomial  of  degree  n, 

<P  = /a2  + t,2.  Spherical  distance  between  P and  Q, 


• • 


A(n-l) 


degree  variance  of  degree  n,  where 


A = 425.28  mgal2, 
B = 24, 


= 0.999617, 


Rs  = Radius  of  the  Bjerhammar  sphere, 
rp  and  r^  = Geocentric  radii  to  the  points  P and  Q. 

A3.  POWER  SPECTRA 

1 . The  Power  Spectrum  Along  a Profile 

Assuming  the  statistics  of  the  local  field  are  known, 
the  power  spectrum  for  observations  along  a profile  can  be  deter- 
mined by  the  Fourier  Transform  of  the  geoid  height  covariance  func- 
tion. Thus,  the  power  spectrum  is 

S(u)  = / C(s)e~1u)S  ds.  (A2 


Since 


= coi(us)  - i sin(ws), 


the  power  spectrum  can  be  expressed  as 


/ C(s)  cos (as ) ds  - i / C(s)  sin(us)  ds,  (A2 


and  snct  C(s)  u an  even  function, 


j C's,  ^in'.os)  us  •••  0. 


so  that  the  expression  for  the  spectrum  reduces  to 


S(ui) 


/ C(s)  cos(u)s)  ds. 


(A5) 


Since  C(s)  is  symmetric,  the  continuous  power  spectrum  along  a 
profile  reduces  to 

co 

S(u))  = 2 / C (s ) cos(ws)  ds,  (A6) 

0 

which  can  also  be  expressed  in  discrete  form  as 

n 1 

f(rA.w)  = As[C(0)  + 2 ^ C(qAs)  cns(^^)  (A7) 

q-i  n 

+ C(nAs)  cos(nt)] , 

where  n equally  spaced  values  of  the  covariance  function  are  used, 
r — 0,  1,2,  ...,n, 


Aw  = 


nAs 


The  notation  for  the  various  power  spectra  to  be  used  is  S — U)  for 
the  spectral  value  between  points  on  profiles  i and  j for  frequency  w, 
and  it  is  assumed  tne  along  track  spectrum  is  identical  for  all  pro- 
files in  a local  area  (i  = j),  or 


’ll 


U / 


= $22(<iG) 


S (w ) . 


(A8) 


2 • The  Across  Track  Power  Spectrum 

The  spectrum  S-|0(w)  of  the  covariance  function 
between  points  or  two  adjacent  orofiles  is  also  the  four  :r  transform 
of  the  covariance  function 


C12(w) 


j C^,(T)e  1(^T  ^T’ 

• Oo 
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Again,  since  c12(t)  is  an  even  function  of  t, 

oo 

S,?(u>)  = 2 / C,?(t)  cos  ( ujt ) d x . 


(A10) 


From  Figure  A1 , if  a is  the  separation  between  the  profiles,  then 

C12(t)  = C(/^T^),  (All) 

so  that  the  crosstrack  spectrum  is 

oo 

S-|2(<o)  = 2 ^ C(  /a2  + t2 ) cos(u>t)c1t.  (A12) 

Profile  1 


Similarly,  the  power  spectrum  Sp^(o))  between  the  point  P and  profile  1, 
Figure  A2  is 


Sp,  (u>)  = 2 / C(  / a?  + t2)  cos(a)t)dT. 

m Q i 

The  discrete  form  for  equation  (A12)  is 

S(rAw)  = As[C1?(0)  + 2 l C^CqAx)  cos 


(A13) 


q=l 


(A14) 


+ C-j^CnAt)  cos  (nr)]. 
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B 

u 


where 

C12(qA-r)  = C(  /a|  + (qAt)2) , 

r * 0,  1,  2,  ....  n, 
q = 0,  1 , 2,  ...»  n. 

The  transformation  of  the  data  from  the  distance  domain  to  the  fre- 
quency domain  is  followed  by  a smoothing  process  to  reduce  spectral 
leakage  and  aliasing  [8]. 

The  power  spectrum  computed  from  equations  (A7)  and  (A14) 
is  plotted  in  Figure  A3.  The  power  spectrum  along  the  profile  is 
represented  by  the  section  in  the  x,z-plane.  The  power  spectrum  S13 
across  the  profile,  from  profile  1 to  profile  3,  is  represented  by 
the  plot  in  the  y,z-plane.  The  power  spectrum  Sp^  of  a point  P is 
represented  by  the  plot  in  a plane  parallel  to  the  x,z-plane  at  some 
distance  along  the  y-axis. 

A4.  GEOID  HEIGHT  PREDICTION 

1 . The  Fundamental  Equation  for  Prediction 

The  predicted  geoid  height  yp  is  assumed  to  be  a linear 
function  of  the  measured  values  x^  so  that  the  fundamental  equation 
for  prediction  [3]  along  one  profile  is 

yp(x)  = / hp(x  - a)  x(a)  da.  (A15) 

— 00 

The  coefficients  hp  depend  on  the  displacements  of  the  measurements 
relative  to  point  P,  not  their  magnitudes.  The  prediction  of  the 
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Finure  A3.  Power  Spectrum  of  Geoid  Heights. 


I 


geoid  height  at  a point  P from  geold  height  data  along  n parallel 
profiles  is  given  by  the  integral  equation 


yp(t)  = l / hpi(x-a)x^ (a)da,  (Alf 

i =1  -® 

where  x^(o)  is  the  geoid  height  measurement  on  the  i—  profile  at 
position  (time)  a,  and  hpi>  (x-a)  is  the  weighting  function  relating 
x.j (a)  to  the  prediction  point  P at  time  x (Figure  A4).  The  contri- 
bution of  n profiles  is  contained  in  the  summation. 


x^a) 


t = a 


t = x 


Profile  i 


Profile  j 


Figure  A4.  Prediction  of  Point  P at  Time  x. 


2.  The  Autocovariance  of  Measurements 


The  measurements  xi  consist  of  measurement  error  (noise) 
n.j  and  true  value  (signal)  s.  so  that 


xi  ■ si  + v 


(A17) 


and  n^  and  si  are  assumed  to  be  uncorrelated  [3].  Then  the  auto- 


covariance function  for  measurements  x on  profiles  i and  j is 

C. . = C . • + D . . , 
ij  U iJ’ 


(A18) 
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where 

C . = the  autocovariance  of  observations  x. , x., 

• \J  ' 

C.  . = the  autocovariance  of  the  signal  s. , s., 

• J J 

D. .  = the  autocovariance  of  noise  n.,  n.. 

• 3 l J 

3.  System  Functions  for  Optimal  Prediction 

In  this  section,  the  Wiener-Hopf  equation  is  solved  for 
the  system  functions  whose  inverse  Fourier  transforms  are  the 
weighting  functions  for  an  optimal  prediction  by  equation  (A16). 

The  geometry  for  prediction  from  two  profiles  is  shown  in  Figure  A 5. 


xi  (t) 


•t  = const. 


«j<'> 


Figure  A5.  Combined  Filtering  and  Prediction. 


The  Wiener-Hopf  system  of  linear  equations  is 

l / hn.(u)  C.  .(u-v)du  = Cp  -(v),  j = 1,  2, 

i=l  -co  : 1 'J 

or  in  the  compact  convolution  notation 


l hp.  (v)*  C (v)  - f.  (v),  .1  = 1,2,  ...  n. 

j =]  r cJ  ' .1 
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Its  spectral  representation  is 
n 


(A20) 


Hpi-  (to)  S..(<o)  - Spj(w),  j = 1,  2,  ...  n 

which  can  easily  be  solved  for  the  system  functions  Hp  . (<o),  provided 
the  spectrums  of  the  covariance  functions  are  known  (Section  A5).  A 
sample  weighting  function  in  the  frequency  domain  as  a function  of 
distance  from  profile  1 is  shown  in  Figure  A6. 

A5.  SYSTEM  FUNCTIONS  FOR  THE  TWO  PROFILE  CASE 

For  prediction  at  point  P between  two  profiles  (Figure  A5), 
the  Wiener-Hopf  equation  (A20)  in  the  frequency  domain  expands  to 


Based  on  the  assumption  that  the  geoid  field  is  stationary  and  the 
profiles  are  parallel  tracks,  it  follows  that  ^ i ( to ) = S^U)  and 

SpU)  = (“)  so  that  the  solution  to  this  system  becomes: 

^1 1 (cj)  Sp-J  (uj)  - S^U)  Sp2(w)^ 


Hpi  (to)  = 


and 


where 


Hp2(w) 


Hd.(u) 

Pi 


S-|  -j  (to) 


S1i(o.)2  - Sj2_(u>)2 


S- 1 (o>)  Sp2  (to ) - S^to)  Sp^(w) 

Si  i (to)2  - S^ 2 (tJ°) 2 


(A23) 


(A24) 


= the  system  weighting  function  for  point  P from 
profile  i , 

= the  qeoid  height  power  spectrum  along  profile, 
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Sp-J  (a;)  = Hp j (to)  Sj  j (to)  + Hp2(to)  S_21  (w) , 

(A21) 

• 

• 

Sp2(to)  = Hp  -j  (to)  S^  2 (to)  + Hp2(<o)  ^22  (to) . 

(A22) 

• • 


• • 


• • 


• • 
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S12(o>)  = the  geoid  height  power  spectrum  across  profile, 

Sp.(oj)  = the  geoid  height  power  spectrum  relating  point  P to 
profile  j. 

A6.  ERROR  ANALYSIS 

^ • The  Error  Covariance  Function 

The  error  of  prediction  Ep  is  the  deviation  of  the 
predicted  value  yp  from  the  true  value  sp, 

eP  = SP  " V (A25) 

or  in  the  continuous  case  (Figure  A7), 


= sp(x)  - y (T). 


Then  by  Equation  (A16) 


(A26) 


n 


sp(T)  ~ l f hp . (t-a)x. (a)da. 


1 = 1 - 


(A  27) 


n °o 


eq(°)  = Sg(o)  - J / hpj (-0)xj(e)d3. 


(A28) 


x.(a) 


ith  Profile 


jth  Profile 


Figure  A7.  Error  Analysis  Between  Two  Points 
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The  statistical  mean,  M,  of  the  product  epeq  is  the  error  covariance 


of  prediction  for  points  P and  Q,  and  is  denoted  by  opq. 


Then, 


°PQ  = M{ep(T)eQ(°)}- 


°pq(T)  = cpn(T)  “ I hPi(T_a)  Cni.  (-a)da 


'PQ 


Pi 


'Qi 


(A29) 


(A30) 


/ hQ.(-B)  Cpj- (r-6)d6 


+ / / hpi-  (i-a)  hg  . ( - e ) .(a-g)dadB, 


-00  -CO 


which  is  best  solved  by  means  of  spectral  representation  [9].  The 
Fourier  transform  of  equation  (A30)  is 


Epq(w)  = spq(w)  - Hpi(u)  Sg.j(w)  - spj(w) 

+ Hpi(u)  %(u)  SjjM- 


(A31 ) 


After  the  weighting  functions  have  been  determined  by  the  procedure 
outlined  in  Section  A5,  the  error  spectrum  is  easily  computed  with 
equation  (A31).  Then  the  error  covariance  function  apq(i)  is  the 
inverse  Fourier  transform  of  Epg;  that  is. 


PQ' 


(t)  2tt  f ^Dr>(u)e 


lWT 


PQ’ 


(A32) 


for  different  combinations  of  P and  Q. 

2 . Error  Variances  for  1°  x 1°  Blocks 

Assuming  a sufficient  distribution  of  predicted  points, 
the  variance  of  the  geoid  height  for  a 1°  x 1°  block  is  the  mean 
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r 


error  covariance  over  all  pairs  of  points  Pi  and  P^  in  the  block; 
that  is. 


, n n 

°2  = ~ I I °i  y 

n2  i=l  j=l  1J 


(A33) 


As  an  example  with  n = 9, 

o-| 2 is  the  error  covariance  between  point  1 (P) 
and  2 (Q) 

073  is  the  error  covariance  between  point  7 (P) 

and  point  3 (Q). 

3 . The  Error  Covariance  Between  1°  x 1°  Blocks 

The  error  covariance  between  1°  x 1°  blocks  I and  II  is 

the  mean  covariance  between  all  pairs  of  points  P.  in  I and  0.  in  II; 

J 

I II 


1 

2 

3 

10 

11 

12 

4 

5 

6 

13 

14 

15 

7 

8 

9 

16 

17 

18 

that  is,  for  n = 9 


1 


9 18 


o(I,  II)  - ^ l l o 
n2  i=l  j=10  1J 


(A34) 


A • The  Error  Integral 

The  error  inteqral,  often  useful  in  studies  of  trafec- 
tory  error  or  gravity  survey  requirements,  is  defined  a:. 


~ j J o(<p,\  <(> ' ,x ' ) 0054)'  d<f.'  dx', 

R2  A 


(A35) 
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m 


I 


h 

where  A is  the  region  of  correlation,  and  <f>‘,  x'  vary  over  A.  For 
the  discrete  case  using  the  error  covariances  defined  above,  this 
becomes 

L(<t>j  ) n 

— — = l o(I,  II)  C0S()>TT  A(Jjtt  AXtt.  (A36) 

r2  ii=i  11  11  11 

A7.  COMPUTATIONAL  RESULTS 

Prediction  error  variances  were  computed  for  profile 
spacings  of  1°,  2°,  3°20'  and  5°.  Geoid  height  sigmas  (a^)  of  0.0, 

0.3,  0.6,  1.2  and  3.0  meters  were  used  in  each  case.  Error  variances 
and  error  integrals  in  meters2  for  the  1°  x 1°  squares  are  shown  in 
Tables  A1  and  A2,  respectively.  The  error  of  prediction  reaches  its 
maximum  value  if  the  computation  point  is  midway  between  two  profiles. 
If  - 3.0  meters  and  the  two  profiles  are  spaced  1°,  2°,  3°20'  or 
5°  apart,  then  the  maximum  errors  of  prediction  are  3.61,  3.85,  4.29 
or  4.89  meters,  respectively. 

The  computed  error  variances  for  prediction  of  geoid  heights 
ai  points  between  profiles  spaced  1°,  2C,  3°20‘  and  5°  apart  are 
plotted  in  Figure  A8.  The  total  error  variance  for  continuous 
prediction  between  profiles  is  the  total  area  under  the  curve.  The 
total  error  variance  increases  as  the  distance  between  profiles 
increases,  although  error  variance  for  prediction  near  a profile 
decreases  <is  Jie  d stance  between  input  profiles  increases. 
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Table  A1 

Prediction  Error  Variances  for  1°  x 1°  Squares 
(meters2) 


Profile 

Spacing 


Table  A2 


Error  Integrals  for  1°  x 1°  Squares 
(meters2) 


Profile 

Spacing 


Geoid  Height  Errors  (one  sigma) 


0.0  0.3  0.6  1.2  3.0 


.011 

.011 

.011 

.011 

.014 

.021 

.021 

.021 

.022 

.022 

.016 

.016 

.016 

.017 

.020 

.013 

.013 

.014 

.014 

.017 

.037 

.037 

.037 

.039 

.048 

.048 

.048 

.048 

.050 

.062 

.037 

.037 

.038 

.039 

.048 

.014 

.014 

.014 

.014 

.017 
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